spatial12
Meredith Franklin
11/22/2019
Packages
library(spacetime) # for s-t analysis
library(gstat) # for s-t analysis (krigST)
library(mgcv)
library(mgcViz)
library(leaflet)
library(sp)
library(RColorBrewer)
library(lubridate)
load("/Users/mf/Dropbox (University of Southern California)/Courses/PM569/Fall 2019/Lecture12/STObj3.rda")
STObj4 <- STObj3[, "1993-07-01::1993-07-31"] #subset to July 1993
STObj4[1,] # the data (date, julian date, year, month, day, id, temp, timeindex)## julian year month day id z timeIndex
## 1993-07-01 728111 1993 7 1 3804 82 1278
## 1993-07-02 728112 1993 7 2 3804 84 1279
## 1993-07-03 728113 1993 7 3 3804 88 1280
## 1993-07-04 728114 1993 7 4 3804 90 1281
## 1993-07-05 728115 1993 7 5 3804 92 1282
## 1993-07-06 728116 1993 7 6 3804 91 1283
## 1993-07-07 728117 1993 7 7 3804 92 1284
## 1993-07-08 728118 1993 7 8 3804 94 1285
## 1993-07-09 728119 1993 7 9 3804 96 1286
## 1993-07-10 728120 1993 7 10 3804 93 1287
## 1993-07-11 728121 1993 7 11 3804 94 1288
## 1993-07-12 728122 1993 7 12 3804 85 1289
## 1993-07-13 728123 1993 7 13 3804 92 1290
## 1993-07-14 728124 1993 7 14 3804 90 1291
## 1993-07-15 728125 1993 7 15 3804 80 1292
## 1993-07-16 728126 1993 7 16 3804 85 1293
## 1993-07-17 728127 1993 7 17 3804 87 1294
## 1993-07-18 728128 1993 7 18 3804 90 1295
## 1993-07-19 728129 1993 7 19 3804 87 1296
## 1993-07-20 728130 1993 7 20 3804 86 1297
## 1993-07-21 728131 1993 7 21 3804 82 1298
## 1993-07-22 728132 1993 7 22 3804 81 1299
## 1993-07-23 728133 1993 7 23 3804 82 1300
## 1993-07-24 728134 1993 7 24 3804 87 1301
## 1993-07-25 728135 1993 7 25 3804 92 1302
## 1993-07-26 728136 1993 7 26 3804 94 1303
## 1993-07-27 728137 1993 7 27 3804 92 1304
## 1993-07-28 728138 1993 7 28 3804 95 1305
## 1993-07-29 728139 1993 7 29 3804 86 1306
## 1993-07-30 728140 1993 7 30 3804 76 1307
## 1993-07-31 728141 1993 7 31 3804 85 1308pm25 <-read.csv("/Users/mf/Dropbox (University of Southern California)/Courses/PM569/Fall 2019/data/pm25_ca_2018.csv")
pm25_sites<-read.csv("/Users/mf/Dropbox (University of Southern California)/Courses/PM569/Fall 2019/data/pm25_sites.csv")
pm<-merge(pm25, pm25_sites, by='site_id')
pm_pal = colorNumeric(c('darkgreen','goldenrod','brown','brown','brown'),
domain=pm$pm25)
leaflet(pm) %>%
addProviderTiles('CartoDB.Positron') %>%
addCircles(label=~paste0(pm25, ' ug/m3'), color=~pm_pal(pm25),
opacity=1, fillOpacity=1, radius=500) %>%
addLegend('bottomleft', pal=pm_pal, values=pm$pm25,
title='PM2.5 (ug/m3)', opacity=1)SpaceTime Variograms
- We use temperature data and construct separable and joint semivariograms
vv <- variogram(object = z ~ lon + lat,
data = STObj4, # July data
width = 80, # spatial bin (80 km)
cutoff = 1000, # consider pts < 1000 km apart
tlags = 0.01:6.01) # 0 days to 6 days
# st separable variogram
sepVgm <- vgmST(stModel = "separable",
space = vgm(10, "Exp", 400, nugget = 0.1),
time = vgm(10, "Exp", 1, nugget = 0.1),
sill = 20)
sepVgm <- fit.StVariogram(vv, sepVgm)
plot(vv,sepVgm,wireframe=TRUE,all=TRUE)
# st joint variogram (metric, see notes)
metricVgm <- vgmST(stModel = "metric",
joint = vgm(10, "Exp", 400, nugget = 0.1),
sill = 20,
stAni = 100)
metricVgm <- fit.StVariogram(vv, metricVgm)
plot(vv,metricVgm,wireframe=TRUE,all=TRUE)
# st joint variogram (sum metric, see notes)
summetricVgm <- vgmST(stModel = "sumMetric",
space = vgm(10, "Exp", 400, nugget = 0.1),
time = vgm(10, "Exp", 1, nugget = 0.1),
joint = vgm(10, "Exp", 400, nugget = 0.1),
sill = 20,
stAni = 100)
summetricVgm <- fit.StVariogram(vv, summetricVgm)
plot(vv,summetricVgm,wireframe=TRUE,all=TRUE)
SpaceTime Kriging
spat_pred_grid <- expand.grid(
lon = seq(-100, -80, length = 20),
lat = seq(32, 46, length = 20)) %>%
SpatialPoints(proj4string = CRS(proj4string(STObj3)))
gridded(spat_pred_grid) <- TRUE
## ------------------------------------------------------------------------
temp_pred_grid <- as.Date("1993-07-01") + seq(3, 28, length = 6)
## ------------------------------------------------------------------------
DE_pred <- STF(sp = spat_pred_grid, # spatial part
time = temp_pred_grid) # temporal part
STObj5 <- as(STObj4[, -14], "STIDF") # convert to STIDF
STObj5 <- subset(STObj5, !is.na(STObj5$z)) # remove missing data
pred_kriged <- krigeST(z ~ lon + lat, # latitude trend
data = STObj5, # data set w/o 14 July
newdata = DE_pred, # prediction grid
modelList = metricVgm, # semivariogram
computeVar = TRUE) # compute variances
color_pal <- rev(colorRampPalette(brewer.pal(11, "Spectral"))(16))
stplot(pred_kriged,
main = "Predictions (degrees Fahrenheit)",
layout = c(3, 2),
col.regions = color_pal)
pred_kriged$se <- sqrt(pred_kriged$var1.var)
stplot(pred_kriged[, , "se"],
main = "Prediction std. errors (degrees Fahrenheit)",
layout = c(3, 2),
col.regions = color_pal)
## ST GAM
- We show a separable model using smooths for space and time as well as a space-time interaction model
- Space is usually defined by a thin plate spline and time as a cubic regression spline. These can be kept as separate components or can be multiplied to create an interactive s-t surface.
- need the mgcv package
pm$date2<-as.Date(pm$date,"%m/%d/%y")
pm$month <- month(pm$date2)
pm$dow<-wday(pm$date2)
pm$jday = as.numeric(format(pm$date2, '%j'))
# separate space and time
gam1 <- gam(pm25~s(longitude,latitude, k=75,bs="tp") + s(jday,k=15,bs="cr") + s(month,k=8,bs="cc"), data=pm, method="ML")
summary(gam1)##
## Family: gaussian
## Link function: identity
##
## Formula:
## pm25 ~ s(longitude, latitude, k = 75, bs = "tp") + s(jday, k = 15,
## bs = "cr") + s(month, k = 8, bs = "cc")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.66460 0.07318 145.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(longitude,latitude) 65.054 71.36 27.34 <2e-16 ***
## s(jday) 13.848 14.00 88.72 <2e-16 ***
## s(month) 5.783 6.00 21.21 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.223 Deviance explained = 22.6%
## -ML = 62973 Scale est. = 91.473 n = 17079plot(gam1)
vgam1<-getViz(gam1)
plot(vgam1)## Hit <Return> to see next plot:
## Hit <Return> to see next plot:
## Hit <Return> to see next plot:
gam2 <- gam(pm25~t2(longitude,latitude,jday, k=c(50,6),d=c(2,1),bs=c("tp","cr")) + s(month, bs="cc",k=8), data=pm, method="ML")
pen.edf(gam2) ## t2(longitude,latitude,jday)rr t2(longitude,latitude,jday)nr
## 125.058140 10.097025
## t2(longitude,latitude,jday)rn s(month)
## 80.156901 5.887925# 3D plot
vis.gam(gam2)
plot(gam2)
gam3 <- gam(pm25~te(longitude,latitude,jday, k=c(10,6),d=c(2,1),bs=c("tp","cr")), data=pm, method="ML")
pen.edf(gam3) ## te(longitude,latitude,jday)1 te(longitude,latitude,jday)2
## 56.93532 56.93532# 3D plot
vis.gam(gam3)